PVA of the NBI: Analysis of Stochastic Event Scenarios - Plots
- Hypothesis: We hypothesize that at the present time the northern bald ibis population can survive without further management and release. We predict that the observed demographic rates will ensure population growth and do not differ between the breeding colonies.
- Study area: Austria, Germany and Italy
- Data: Results of the second NetLogo model and tables created in script 07b
In this script we create barplots to show the distribution of survival and reproduction values between the scenarios per combination of stochastic event frequency and severity and we show the distribution between the number of supplements and the years of supplementation for scenarios with positive population growth.
Setup
## for non-CRAN packages please keep install instruction
## but commented so it is not run each time, e.g.
# devtools::install_github("EcoDynIZW/template")
## libraries used in this script
## please add ALL LIBRARIES NEEDED HERE
## please remove libraries from the list that are not needed anymore
## at a later stage
library("here")
library("tidyverse")
library("magrittr")
library("ggplot2")
library("viridisLite")
library("gridExtra")##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
##
## combine
Data
Barplot of all scenarios with positive population growth
Version 1: Absolute numbers
Count the individuals a) Overall
##
## 0.2 0.3 0.36
## 160 160 500
##
## 0.08 0.19 0.26
## 160 160 500
##
## 0.14 0.24 0.31
## 160 160 500
##
## 0.02 0.14 0.22
## 240 240 340
##
## 0.53 0.58 0.66 1.06 1.41 3.97
## 400 80 80 80 80 100
- Scenarios with positive population growth
##
## 0.2 0.3 0.36
## 160 160 394
##
## 0.08 0.19 0.26
## 160 160 394
##
## 0.14 0.24 0.31
## 160 160 394
##
## 0.02 0.14 0.22
## 240 214 260
##
## 0.53 0.58 0.66 1.06 1.41 3.97
## 294 80 80 80 80 100
Create a matrix
StevS_mort_RR_mat_all <- matrix(data = c(
500,160,160,NA, # Baseline, 10%, 25%
500,160,160,NA, # Baseline, 10%, 25%
500,160,160,NA, # Baseline, 10%, 25%
340,240,240, NA, # Baseline, 10%, 25%
400,80,80,80,80,100 # Baseline, 10%, 25%, Status quo, All chicks
), byrow = F)
StevS_mort_RR_mat_all## [,1]
## [1,] 500
## [2,] 160
## [3,] 160
## [4,] NA
## [5,] 500
## [6,] 160
## [7,] 160
## [8,] NA
## [9,] 500
## [10,] 160
## [11,] 160
## [12,] NA
## [13,] 340
## [14,] 240
## [15,] 240
## [16,] NA
## [17,] 400
## [18,] 80
## [19,] 80
## [20,] 80
## [21,] 80
## [22,] 100
StevS_mort_RR_mat <- matrix( data = c(
394,160,160,NA,
394,160,160,NA,
394,160,160,NA,
260,214,240,NA,
294,80,80,80,80,100
), byrow = F)
StevS_mort_RR_mat## [,1]
## [1,] 394
## [2,] 160
## [3,] 160
## [4,] NA
## [5,] 394
## [6,] 160
## [7,] 160
## [8,] NA
## [9,] 394
## [10,] 160
## [11,] 160
## [12,] NA
## [13,] 260
## [14,] 214
## [15,] 240
## [16,] NA
## [17,] 294
## [18,] 80
## [19,] 80
## [20,] 80
## [21,] 80
## [22,] 100
Plot
par(mar = c(5,6,1,1))
barplot(StevS_mort_RR_mat_all, beside = T, ylim = c(0,700), space = c(0,-0.3), las = 1,
col = "grey", cex.axis = 1.5, cex.lab = 1.5, border = NA)
barplot(StevS_mort_RR_mat, beside = T, ylim = c(0,700), space = c(0,-0.3), las = 1,
col = c(viridis(4, end = 0.95), viridis(4, end = 0.95), viridis(4, end = 0.95), viridis(4, end = 0.95), viridis(4, end = 0.95), cividis(2, begin = 0.3, end = 0.9)), cex.axis = 1.5, cex.lab = 1.5, ylab = "", border = NA, add = T)
axis(1, at = c(1.2, 4.9,8.6,12.3,16.5), labels = c("","","","",""), cex.axis = 1.5)
axis(1, at = c(1.2, 4.9,8.6,12.3,16.5), labels = c("s1","s2","s3","s4","RR"), cex.axis = 1.5, line = 0.2, lwd = 0)
mtext("Count", side = 2, line = 4, cex = 1.5)
box()Version 2: Percentage
Calculate the percent of the used s1-s4 and reproduction values in scenarios with positive population growth
## [,1]
## [1,] 78.80000
## [2,] 100.00000
## [3,] 100.00000
## [4,] NA
## [5,] 78.80000
## [6,] 100.00000
## [7,] 100.00000
## [8,] NA
## [9,] 78.80000
## [10,] 100.00000
## [11,] 100.00000
## [12,] NA
## [13,] 76.47059
## [14,] 89.16667
## [15,] 100.00000
## [16,] NA
## [17,] 73.50000
## [18,] 100.00000
## [19,] 100.00000
## [20,] 100.00000
## [21,] 100.00000
## [22,] 100.00000
Plot
par(mar = c(4,8,1,1))
barplot(StevS_mort_RR_mat_perc, beside = T, ylim = c(0,100), space = c(0,-0.3), las = 1,
col = c(viridis(4, end = 0.95), viridis(4, end = 0.95), viridis(4, end = 0.95), viridis(4, end = 0.95), viridis(4, end = 0.95), cividis(2, begin = 0.3, end = 0.9)), cex.axis = 2.5, cex.lab = 3, ylab = "", border = NA)
axis(1, at = c(1.2, 5.1,9.2,13.2,18.7), labels = c("","","","",""), cex.axis = 2.5)
axis(1, at = c(1.2, 5.1,9.2,13.2,18.7), labels = c("s1","s2","s3","s4","RR"), cex.axis = 2.5, line = 0.2, lwd = 0)
mtext("% of used values in scenarios", side = 2, line = 6, cex = 3)
mtext("with positive population growth", side = 2, line = 4, cex = 3)
box()
text(x = -2.8, y = -5, labels = "b)", xpd=NA, cex = 3)Save plot
# png(filename = here("plots", "04_PVA", "Stoch_eve_Supplements_bar.png"), width = 1200, height = 772)
#
# par(mar = c(4,8,1,1))
# barplot(StevS_mort_RR_mat_perc, beside = T, ylim = c(0,100), space = c(0,-0.3), las = 1,
# col = c(viridis(4, end = 0.95), viridis(4, end = 0.95), viridis(4, end = 0.95), viridis(4, end = 0.95), viridis(4, end = 0.95), cividis(2, begin = 0.3, end = 0.9)), cex.axis = 2.5, cex.lab = 3, ylab = "", border = NA)
# axis(1, at = c(1.2, 5.1,9.2,13.2,18.7), labels = c("","","","",""), cex.axis = 2.5)
# axis(1, at = c(1.2, 5.1,9.2,13.2,18.7), labels = c("s1","s2","s3","s4","RR"), cex.axis = 2.5, line = 0.2, lwd = 0)
# mtext("% of used values in scenarios", side = 2, line = 6, cex = 3)
# mtext("with positive population growth", side = 2, line = 4, cex = 3)
# box()
# text(x = -2.8, y = -5, labels = "b)", xpd=NA, cex = 3)
#
# dev.off()A combinded figure of the barplots for the Baseline and Improvement scenarios and the stochastic event and supplement scenarios
# png(filename = here("plots", "04_PVA", "combined_bar.png"), width = 1200, height = 772)
# par(mar = c(0.4,4.6,1,0) + 0.1,
# mfrow= c(2,1),
# oma = c(3, 4, 0, 0) + 0.1)
# #mai = c(1, 0.9, 0, 0.1))
#
# barplot(mort_RR_mat_perc, beside = T, ylim = c(0,100), space = c(0,-0.3), las = 1,
# col = viridis(4, end = 0.95), cex.axis = 2.5, cex.lab = 3, border = NA,
# ylab = " ")
# mtext("% of used values in scenarios", side = 2, line = 6.7, cex = 3, adj = 2)
# mtext("with positive population growth", side = 2, line = 4.7, cex = 3, adj = 2)
# box()
# text(x = -0.7, y = 90, labels = "a)", xpd=NA, cex = 3)
# #text(x = -2.8, y = -5, labels = "positive", xpd=NA, cex = 3)
#
# barplot(StevS_mort_RR_mat_perc, beside = T, ylim = c(0,100), space = c(0,-0.3), las = 1,
# col = c(viridis(4, end = 0.95), viridis(4, end = 0.95), viridis(4, end = 0.95), viridis(4, end = 0.95), viridis(4, end = 0.95), cividis(2, begin = 0.3, end = 0.9)), cex.axis = 2.5, cex.lab = 3, ylab = "", border = NA)
# axis(1, at = c(1.2, 5.1,9.2,13.2,18.7), labels = c("","","","",""), cex.axis = 2.5, outer = F)
# axis(1, at = c(1.2, 5.1,9.2,13.2,18.7), labels = c("s1","s2","s3","s4","RR"), cex.axis = 2.5, line = 0.2, lwd = 0, outer = F)
# #mtext("% of used values in scenarios", side = 2, line = 6, cex = 3)
# #mtext("with positive population growth", side = 2, line = 4, cex = 3)
# box()
# text(x = -0.7, y = 90, labels = "b)", xpd=NA, cex = 3)
# dev.off()Stochastic event frequency and severity
df_StevS_summ <- df_StevS_lambda_pos %>%
group_by(Stoch_Frequency, Stoch_Severity) %>%
summarise(Freq = n())##
## 0.05 0.1 0.15 0.2 0.25
## 0.05 41 41 41 41 41
## 0.1 41 41 41 41 41
## 0.15 41 41 41 41 41
## 0.2 41 41 41 41 41
##
## 33 34 36 37
## 5 1 3 11
## [1] 714
df_StevS_summ$Freq_perc <- (df_StevS_summ$Freq/41)*100 # value of the number of scenarios with positive population growth per combination of stochastic event severity and frequency / the overall number of scenarios per combination of stochastic event frequency and severity * 100 for percentage
range(df_StevS_summ$Freq_perc) # 80.4878 90.2439, the range is not so wide## [1] 80.4878 90.2439
Version 1: Absolute numbers
# Basic Plot
ggplot(data = df_StevS_summ, aes(x = Stoch_Severity,
y = Stoch_Frequency,
fill = Freq)) +
geom_tile() +
# Colour and borders
scale_fill_viridis_c(option = "viridis") +
geom_hline(yintercept = c(0.025, 0.075, 0.125, 0.175, 0.225),
colour = "white",
size = 1) +
geom_vline(xintercept = c(0.025, 0.075, 0.125, 0.175, 0.225, 0.275),
colour = "white",
size = 1) +
# Axis breaks
scale_x_continuous(#name = "Stochastic event severity",
breaks = c(0.05, 0.10, 0.15, 0.20, 0.25), expand = c(0, 0)) +
scale_y_continuous(#name = "Stochastic event frequency",
breaks = c(0.05, 0.10, 0.15, 0.20), expand = c(0, 0)) +
# Theme
theme(rect = element_blank(),
axis.text = element_text(size = 14),
axis.ticks = element_line(size = 1),
axis.ticks.length = unit(.25, "cm"),
axis.title = element_text(size = 16),
legend.title = element_text(size = 16),
legend.text = element_text(size = 14)) +
labs(x = "Stochastic event severity", y = "Stochastic event frequency", fill = "Frequency")Version 2: Percentage
# Basic Plot
ggplot(data = df_StevS_summ, aes(x = Stoch_Severity,
y = Stoch_Frequency,
fill = Freq_perc)) +
geom_tile() +
# Colour and borders
scale_fill_viridis_c(option = "viridis") +
geom_hline(yintercept = c(0.025, 0.075, 0.125, 0.175, 0.225),
colour = "white",
size = 1) +
geom_vline(xintercept = c(0.025, 0.075, 0.125, 0.175, 0.225, 0.275),
colour = "white",
size = 1) +
# Axis breaks
scale_x_continuous(#name = "Stochastic event severity",
breaks = c(0.05, 0.10, 0.15, 0.20, 0.25), expand = c(0, 0)) +
scale_y_continuous(#name = "Stochastic event frequency",
breaks = c(0.05, 0.10, 0.15, 0.20), expand = c(0, 0)) +
# Theme
theme(rect = element_blank(),
axis.text = element_text(size = 14),
axis.ticks = element_line(size = 1),
axis.ticks.length = unit(.25, "cm"),
axis.title = element_text(size = 16),
legend.title = element_text(size = 16),
legend.text = element_text(size = 14)) +
labs(x = "Stochastic event severity", y = "Stochastic event frequency", fill = "Percent") +
geom_text(mapping = aes(y = 0.05, x = 0.05), label = "a)", vjust = 4.8, hjust = 5, size = 6) +
coord_cartesian(clip = "off")Number of supplements & years of supplementation
##
## 0 15 30
## 4 20 200 200
## 7 0 200 200
df_Supp_summ <- df_StevS_lambda_pos %>%
group_by(N_Supplements, Supplement_Time) %>%
summarise(Freq = n())
df_Supp_summ## # A tibble: 5 x 3
## # Groups: N_Supplements [3]
## N_Supplements Supplement_Time Freq
## <dbl> <dbl> <int>
## 1 0 4 20
## 2 15 4 171
## 3 15 7 174
## 4 30 4 174
## 5 30 7 175
Remove the combination with 0 supplements
## # A tibble: 4 x 3
## # Groups: N_Supplements [2]
## N_Supplements Supplement_Time Freq
## <dbl> <dbl> <int>
## 1 15 4 171
## 2 15 7 174
## 3 30 4 174
## 4 30 7 175
##
## 4 7
## 0 20 0
## 15 200 200
## 30 200 200
df_Supp_summ$Freq_perc <- (df_Supp_summ$Freq/200) *100
range(df_Supp_summ$Freq_perc) # 85.5 87.5, very small range## [1] 85.5 87.5
df_Supp_summ$cat <- c(1:4)
cols <- c("N_Supplements", "Supplement_Time", "cat")
df_Supp_summ[cols] <- lapply(df_Supp_summ[cols], as.factor)# Plot with percent of scenarios
ggplot(data = df_Supp_summ, aes(x = Supplement_Time,
y = N_Supplements,
fill = Freq_perc)) +
geom_tile() +
# Colour and borders
scale_fill_viridis_c(option = "viridis") +
geom_hline(yintercept = c(0.5, 1.5, 2.5, 3.5),
colour = "white",
size = 1, na.rm = T) +
geom_vline(xintercept = c(0.5, 1.5, 2.5),
colour = "white",
size = 1) +
# Axis breaks
scale_x_discrete(#name = "Time of supplementation",
breaks = c(4,7), expand = c(0, 0)) +
scale_y_discrete(#name = "Number of supplements",
breaks= c(0,15,30), expand = c(0,0)) +
# Theme
theme(rect = element_blank(),
axis.text = element_text(size = 14),
axis.ticks = element_line(size = 1),
axis.ticks.length = unit(.25, "cm"),
axis.title = element_text(size = 16),
legend.title = element_text(size = 16),
legend.text = element_text(size = 14), legend.box.spacing = unit(2.5, units = "cm")) +
labs(x = "Year of supplementing", y = "Number of supplements", fill = "Percent")+
geom_text(mapping = aes(y = 1, x = 1), label = "b)", vjust = 6.5, hjust = 5.8, size = 6) +
coord_cartesian(clip = "off")Save plot
# ggsave(filename = here("plots", "04_PVA", "Supplements_percent.png"),
# width = 13, height = 8, units = "cm")Session Info
## DO NOT REMOVE!
## We store the settings of your computer and the current versions of the
## packages used to allow for reproducibility
Sys.time()## [1] "2022-01-06 20:15:50 CET"
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22000)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=German_Germany.1252 LC_CTYPE=German_Germany.1252 LC_MONETARY=German_Germany.1252 LC_NUMERIC=C LC_TIME=German_Germany.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] gridExtra_2.3 viridisLite_0.3.0 data.table_1.13.2 magrittr_1.5 forcats_0.5.0 stringr_1.4.0 dplyr_1.0.2 purrr_0.3.4 readr_1.4.0
## [10] tidyr_1.1.2 tibble_3.0.3 ggplot2_3.3.2 tidyverse_1.3.0 here_0.1
##
## loaded via a namespace (and not attached):
## [1] httr_1.4.2 pkgload_1.1.0 jsonlite_1.7.1 showtext_0.9 modelr_0.1.8 assertthat_0.2.1 showtextdb_3.0 blob_1.2.1 cellranger_1.1.0
## [10] yaml_2.2.1 remotes_2.2.0 sessioninfo_1.1.1 pillar_1.4.6 backports_1.1.10 glue_1.4.2 digest_0.6.25 rvest_0.3.6 colorspace_1.4-1
## [19] htmltools_0.5.0 pkgconfig_2.0.3 devtools_2.3.2 broom_0.7.4 haven_2.3.1 bookdown_0.20 sysfonts_0.8.1 scales_1.1.1 processx_3.4.4
## [28] d6_0.1.0.0 farver_2.0.3 generics_0.0.2 usethis_1.6.3 ellipsis_0.3.1 withr_2.3.0 cli_2.0.2 crayon_1.3.4 readxl_1.3.1
## [37] memoise_1.1.0 evaluate_0.14 ps_1.4.0 fs_1.5.0 fansi_0.4.1 xml2_1.3.2 pkgbuild_1.1.0 tools_4.0.3 prettyunits_1.1.1
## [46] hms_0.5.3 lifecycle_0.2.0 munsell_0.5.0 reprex_0.3.0 callr_3.5.0 compiler_4.0.3 tinytex_0.26 rlang_0.4.7 grid_4.0.3
## [55] rstudioapi_0.11 labeling_0.3 rmarkdown_2.4 testthat_2.3.2 gtable_0.3.0 DBI_1.1.0 R6_2.4.1 lubridate_1.7.9 knitr_1.30
## [64] utf8_1.1.4 rprojroot_1.3-2 desc_1.2.0 stringi_1.5.3 rmdformats_0.3.7 Rcpp_1.0.5 vctrs_0.3.4 dbplyr_1.4.4 tidyselect_1.1.0
## [73] xfun_0.18